library(ggplot2)
require(data.table)
library("dplyr")
library("ggpubr") 
require(tidyverse) #for easy data manipulation and visualization
require(caret) #for easy machine learning workflow
library(ncdf4)



##=========directory to replication data and scripts=============
wd1 <- "C:\\Users\\yfan\\...\\NFood_replication\\"
case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")
casename = c("RCP45",expression('RCP45'[(SSP5LU)]),"RCP85",expression('RCP85'[(RCP45CO2)]),"RCP85+SAI","RCP85+MSB","RCP85+CCT")


#******************************************************************************
#=============================Statistics and Figures===========================
#******************************************************************************

#=======load data prepared by Data.Fig1.R=============
load(paste0(wd1,"/data/Rdata/Fig1.GlobalAverage-GrowingSeasonClimate.Rdata")) #5-year running mean


#======global warming by 2050/2100 under different scenarios
#2050(2046-2055) average compared to piControl (1850-1879)
warming2050 <- TSA.gl[Year>=2046 & Year<=2055, mean(tsa.gl, na.rm = T)-piControl.tsa.gl, by=.(Case)]

#2075-2099 average compared to piControl (1850-1879)
warming2090 <- TSA.gl[Year>=2075 & Year<=2099, mean(tsa.gl, na.rm = T)-piControl.tsa.gl, by=.(Case)]
warming2090 #annual global average warming ~5.4 under RCP85

#====growing season warming over cropping area
cropwarming2050= df.5case.wm[Year>=2046 & Year<=2055, mean(tsa.gl)- clim1850.GR.wm$tsa.gl, by=.(Case)]  #>3.5°C warming in RCP85 by 2050

(cropwarming2090= df.5case.wm[Year>=2075 & Year<=2099, mean(tsa.gl)- clim1850.GR.wm$tsa.gl, by=.(Case)]) #crop area warming 6°C under RCP85 
#======compare with IPCC AR5:
#According to IPCC AR5, the global mean warming over land by 2081–2100 is projected to be 3.4~6.2 °C with reference to the 1986-2005 period, 
#which by itself has warmed 0.61 °C since 1850-1900. Thus, total warming over land will be at least 4~6.8 °C above the preindustrial level.
#=======cropwarming2100, 6°C under RCP85, is within the range 4~6.8 °C

#====change in direct and diffuse radiation in SAI 
df.5case.wm[Year>=2075 & Year<=2099, mean(radd.gl), by=.(Case)] #SAI-RCP85: -11.2 W/m2 or -12%
df.5case.wm[Year>=2075 & Year<=2099, mean(radi.gl), by=.(Case)] #SAI-RCP85: +7.2 W/m2 or 20%


#===average cooling in 2050-2069=====
df.5case.wm[Year>=2050 & Year<=2069, .(mean(tsa.gl)), by=.(Case)] #SAI - RCP85:-0.89 degC
df.5case.wm[Year>=2075 & Year<=2099, .(mean(tsa.gl)), by=.(Case)] #SAI - RCP85: -1.8



#save two color-blind-friendly palettes, one with gray, and one with black
#http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/

#install.packages("viridis")
library(viridis)
cbPalette <- viridis(10)[seq(1,10,2)]
barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors
casecolor <- c("#26828EFF", "#440154FF", "#CC6600", "#56B4E9","#B4DE2CFF") #"#35B779FF", 

#casecolor <- c("#252525", "#7E6148FF", "#CC79A7", "#8491B4FF", "#91D1C2FF")
casecolor <- c("#440154FF", "#7E6148FF", "#CC79A7", "#8491B4FF", "#91D1C2FF")
barplot(height=rep(1,length(casecolor)),col=casecolor) #display colors


#============plot global average growing season climate trends over cropping areas
df.clim <- df.5case.ma[, c("Case", "Year", "tsa.gl", "ppt.gl", "radd.gl", "radi.gl", "rh.gl")] #5-year running mean
df.clim$rsds.gl <- df.clim$radd.gl+df.clim$radi.gl #surface total downwelling shortwave radiation
df.clim$Var <- "Var"
df.clim <- df.clim[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Var", by=.(Case,Year)]
names(df.clim)[3:4] <- c("Climate", "Data.ma") 

#==preindustrial baseline climate (growing season mean)
piControl<- data.table(Climate=c("tsa.gl", "ppt.gl","radd.gl", "radi.gl", "rh.gl", "rsds.gl"),
                       Baseline=c(clim1850.GR.wm$tsa.gl, NA,NA,NA,NA,NA))
                       #label=c("Preindustrial baseline", NA,NA,NA,NA,NA))
df.clim <- df.clim[piControl, on=.(Climate)]



##======================Final Separate plot for Temperature, Radiation, Moisture
#show both crop area growing season climate change and background global annual climate change for temperature
str(TSA.gl.ma)  #global annual average climate (without weighting by crop area and growing season)
TSA.gl.ma$Climate <- "tsa.background"
TSA.gl.ma$Baseline <- piControl.tsa.gl  #global annual mean climate of 1850-1879
TSA.gl.ma[,Data.ma:=tsa.gl]; TSA.gl.ma[, tsa.gl:=NULL]


#piControl growing season average climate, assuming same total crop area as in 2006
clim1850.GR.wm$tsa.gl
df.tsa <- df.clim[Climate=="tsa.gl"] 
df.tsa <- rbind(df.tsa, TSA.gl.ma)



cases <- c("RCP85", "RCP45", "SAI", "MSB", "CCT")
override.linetype <- c(1,1,1,1,1)
override.size <- c(0.6, 0.6, 0.5, 0.5,0.5)

##===Temperature===
library(egg) #tag subplots easily but will remove facet strips

facet_text <- data.frame(
    Climate = c("tsa.background", "tsa.gl"), 
    labels = c("Global land annual mean", "Growing season mean"))

p0 <- df.tsa %>% 
  mutate(Case = factor(Case, levels=cases)) %>%
  mutate(Climate = factor(Climate, levels=c("tsa.background", "tsa.gl"))) %>%
  ggplot(aes(Year, Data.ma)) + 
  geom_line(aes(color = Case, linetype = Case, size=Case )) + 
  geom_hline(aes(yintercept=Baseline), linetype="longdash", color="#000000", size = 0.8 ) +
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) +
  scale_colour_manual(name = "Case", labels = cases,  values = casecolor, guide = FALSE) +
  scale_x_continuous(name=NULL, breaks =c(2025, 2050, 2075, 2099)) + #xlab("")+ylab("") +
  scale_y_continuous(name="°C", breaks = scales::pretty_breaks(7)) + #set number of y breaks optimized for each facet  
  facet_wrap(~Climate, ncol=2, dir="v",  #dir controls the vertical or horizontal sequence of facets 
             scales='free_y', strip.position = "left", #use left to edit facet y titles
             labeller= label_parsed) 
#set x = -Inf, y = -Inf, otherwise geom_text will require labels for each Case and Var group
p0 <- p0 + geom_text(data = facet_text, aes(x = -Inf, y = -Inf, hjust = -0.05, vjust = -14.5, label=labels), size=4, color="gray3")
#merge the legends for color and size 
p0 <- p0 + guides(colour = guide_legend(keywidth = 1.5, keyheight=0.8, label.position = "right", label.hjust = 0, ncol=1,
                                        override.aes = list(size = override.size )))
##remove grid and background or remove all
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + 
  ggtitle(expression('('*bold(b)*') Air temperature')) + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=12)) +
  theme(legend.text = element_text(size=11), legend.title = element_blank(), legend.position = c(0.1, 0.63), legend.background = element_blank())+#legend.position = "right") + 
  theme(axis.title=element_text(size=13), axis.text = element_text(size=12)) +
  theme(panel.spacing = unit(0.5, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 11))#, strip.background = element_blank(), strip.placement = "outside")

#tag_facet(p0) #default tag for each subplot
tag_facet(p0, #tag_facet replace facet strips with in-panel tags
          x = -Inf, y = -Inf, vjust = -1.2, hjust = -0.05,
          open = "", close = "", # style for "(" and ")"
          fontface = 3, size = 3.6,  #family = "serif",
          tag_pool = c("preindustrial baseline 6.9 °C",
                       "preindustrial baseline 20.7 °C"))

#or use ggsave to save the last plot on the screen
ggsave(paste0(wd1,"/figures/main/Fig1b.GlobalTrend_TSA.svg"), device ="svg", units = "in", width = 6.5, height = 2.5, dpi = 600) #horizontal panel



##===Radiation===
#use package facetscales to set different limits for different facets
#devtools::install_github("zeehio/facetscales")
#library(g)
library(facetscales)

# scales_y <- list(
#   radd.gl = scale_y_continuous(limits = c(78, 96), breaks = seq(78, 94, 4)),
#   radi.gl = scale_y_continuous(limits = c(34, 44), breaks = seq(34, 42, 2)) )

p0 <-  df.clim[Climate=="radd.gl" | Climate=="radi.gl"] %>%
  #mutate(Climate = factor(Climate, levels=vars, labels=varunit)) %>%
  mutate(Climate = factor(Climate, levels=c("radd.gl", "radi.gl"))) %>%
  mutate(Case = factor(Case, levels=cases)) %>%
  ggplot(aes(Year, Data.ma)) + 
  geom_line(aes(color = Case, size=Case) ) + 
  #geom_hline(aes(yintercept=Baseline), linetype="dotted", color="#000000", size = 0.8 ) +
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) +
  scale_colour_manual(name = "Case", labels = cases,  values = casecolor, guide = FALSE) +
  facet_wrap(~Climate, ncol=2, dir="v",  #dir controls the vertical or horizontal sequence of facets
             scales='free_y', strip.position = "left", labeller= label_parsed) + #facet_wrap does not work with tag_facet
 # facet_grid_sc(rows=vars(Climate), scales = list(y = scales_y)) + #set different scales for subplots (only works for vertical panel)
  scale_x_continuous(name=NULL, breaks =c(2025, 2050, 2075, 2099)) + 
  scale_y_continuous(expand = expand_scale(mult = c(0.05, 0.1), add = c(0, 0))) + #add expansion factors or constants to lower and upper limits
  xlab("") + ylab(expression("W m"^-2*"")) #when scale_y_continuous does not work with facet_grid_sc
  #scale_y_continuous(name=expression("W m"^-2*""), breaks = scales::pretty_breaks(6)) #set number of y breaks optimized for each facet  

#merge the legends for color and size 
#p0 <- p0 + guides(colour = guide_legend(keywidth = 1, label.position = "right", label.hjust = 0, ncol=1, override.aes = list(size = override.size)))
##remove grid and background or remove all
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  ggtitle(expression('('*bold(c)*') Solar radiation')) + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=12)) +
  #theme(legend.text = element_text(size=11), legend.title = element_blank(), legend.position = "right") + 
  theme(axis.title=element_text(size=13), axis.text = element_text(size=12)) +
  theme(panel.spacing = unit(0.5, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 11))#, strip.background = element_blank(), strip.placement = "outside")

tag_facet(p0, #tag_facet replace facet strips with in-panel tags
          x = -Inf, y = -Inf, vjust = -14.3, hjust = -0.05,
          open = "", close = "", # style for "(" and ")"
          fontface = 1, size = 4,  family = "sans",# "serif",
          tag_pool = c("Direct radiation", "Diffuse radiation"))

ggsave(paste0(wd1,"/figures/main/Fig1c.GlobalTrend_Radiation.svg"), device ="svg", units = "in", width = 6.5, height = 2.5, dpi = 600) #horizontal panel


#==surface total downwelling shortwave radiation trend (for R2)
p0 <-  df.clim[Climate=="rsds.gl"] %>%
  mutate(Case = factor(Case, levels=cases)) %>%
  ggplot(aes(Year, Data.ma)) + 
  geom_line(aes(color = Case, size=Case) ) + 
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) +
  scale_colour_manual(name = "Case", labels = cases,  values = casecolor, guide = FALSE) +
  scale_x_continuous(name=NULL, breaks =c(2025, 2050, 2075, 2099)) + 
  scale_y_continuous(expand = expand_scale(mult = c(0.05, 0.1), add = c(0, 0))) + #add expansion factors or constants to lower and upper limits
  xlab("") + ylab(expression("W m"^-2*"")) #when scale_y_continuous does not work with facet_grid_sc
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  ggtitle(expression('('*bold(b)*') Surface downwelling shortwave radiation')) + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=8)) +
  theme(axis.title=element_text(size=9), axis.text = element_text(size=8)) +
  theme(panel.spacing = unit(0.5, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 7))#, strip.background = element_blank(), strip.placement = "outside")
p0 + guides(colour = guide_legend(keywidth = 1.5, keyheight=0.8, label.position = "right", label.hjust = 0, ncol=1,
                                  override.aes = list(size = override.size )))+
  theme(legend.text = element_text(size=8), legend.title = element_blank(), legend.position = c(0.2, 0.75), legend.background = element_blank()) 
  
ggsave(paste0(wd1,"/figures/supplementary/Fig.R2.GlobalTrend_TotRadiation-5yrma.svg"), device ="svg", units = "in", width = 4., height = 2.6, dpi = 600) #horizontal panel



##===Precipitation and RH===
#facet_grid_sc(rows=vars(Climate)..only works for vertical panel 
# scales_y <- list(
#   ppt.gl = scale_y_continuous(limits = c(100, 125), breaks = seq(100, 125, 5)),
#   rh.gl = scale_y_continuous(limits = c(58, 64), breaks = seq(58, 63, 1)) )
#==for horizontal panel use expand_scale to adjust the y limits to zoom

p0 <-  df.clim[Climate=="ppt.gl" | Climate=="rh.gl"] %>% 
  mutate(Climate = factor(Climate, levels=c("ppt.gl", "rh.gl"))) %>%
  mutate(Case = factor(Case, levels=cases)) %>%
  ggplot(aes(Year, Data.ma)) + 
  geom_line(aes(color = Case, size=Case) ) + 
  #geom_hline(aes(yintercept=Baseline), linetype="dotted", color="#000000", size = 0.8 ) +
  scale_linetype_manual(values=override.linetype, guide = FALSE)+ 
  scale_size_manual(values=override.size, guide = FALSE) +
  scale_colour_manual(name = "Case", labels = cases,  values = casecolor, guide = FALSE) +
  facet_wrap(~Climate, ncol=2, dir="v",
             scales='free_y', strip.position = "left", labeller= label_parsed) +
  #facet_grid_sc(rows=vars(Climate), scales = list(y = scales_y)) +
  scale_x_continuous(name=NULL, breaks =c(2025, 2050, 2075, 2099)) + 
  scale_y_continuous(expand = expand_scale(mult = c(0.05, 0.15), add = c(0, 0)))+  #add expansion factors or constants to lower and upper limits 
  #xlab("")+ ylab(expression("         %                             mm month"^-1*""))
  ylab(expression('mm month'^-1))
   
#merge the legends for color and size 
# p0 <- p0 + guides(colour = guide_legend(keywidth = 1, label.position = "right", label.hjust = 0, ncol=1, override.aes = list(size = override.size)))
##remove grid and background or remove all
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  ggtitle(expression('('*bold(d)*') Moisture')) + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=12)) +
  #theme(legend.text = element_text(size=13), legend.title = element_blank(), legend.position = "right") + 
  theme(axis.title=element_text(size=13), axis.text = element_text(size=12)) +
  theme(panel.spacing = unit(0.5, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 11))#, strip.background = element_blank(), strip.placement = "outside")

p0 <- tag_facet(p0, #tag_facet replace facet strips with in-panel tags
          x = -Inf, y = -Inf, vjust = c(-11., -10.8), hjust = -0.07,
          open = "", close = "", # style for "(" and ")"
          fontface = 1, size = 4,  family = "sans",# "serif", 
          tag_pool = c("'Precipitation'", "'Relative Humidity (%)'"), parse = TRUE)
          #tag_pool = c(expression("'Precipitation (mm month'^-1*')'"), "'Relative Humidity (%)'"), parse = TRUE)
p0
ggsave(paste0(wd1,"/figures/main/Fig1d.GlobalTrend_Moisture.svg"), device ="svg", units = "in", width = 6.7, height = 2.5, dpi = 600) #horizontal panel





#=======check synchronous drop in PPT and RH and peak of RH effect around year 2067=====
df.5case.wm[Case=="RCP85" & Year>=2065 & Year<=2075, .(Year, ppt.gl, rh.gl), by=.(Case)] 

par(mar=c(1, 5, 1, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
split.screen(c(2,1))
screen(1) # prepare screen 1 for output
plot(ppt.gl~Year, data=df.5case.wm[Case=="RCP85"], type="l", bty="n", xaxt="n", xlab="", ylab="Precipitation (mm/month)")
screen(2)
par(mar=c(2.3, 5, 1, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(rh.gl~Year, data=df.5case.wm[Case=="RCP85"], type="l", bty="n", xaxt="n", xlab="", ylab="Relative humidity (%)")
close.screen(all = TRUE)

par(mar=c(2.3, 4.6, 0.5, 1), mgp=c(1.4,0.5,0), las=1, new=TRUE) #overlay plot to the above split plot
plot(ppt.gl~Year, data=df.5case.wm[Case=="RCP85"], type="n", bty="n", xaxt="l", yaxt="n", 
     at=c(2025, 2050, 2075, 2099), xlab="Year", ylab="", cex=0.9)
abline(v=2066, lty="dotted") 
text(x=2066, y=110, "2066", pos=2, srt=90) #mark the real year when synchronous drop in P and RH happens
text(x=2020, y=125, "RCP85", pos=2, srt=0) 


#=======add label for year 2067 for responses letter======
p0 <- p0 + geom_vline(aes(xintercept=2067), linetype="dotted", color="#000000", size = 0.8 ) 
p0 <- p0 + guides(colour = guide_legend(keywidth = 1, label.position = "right", label.hjust = 0, ncol=1, override.aes = list(size = override.size)))
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  ggtitle(expression('('*bold(d)*') Moisture')) + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=12)) +
  theme(legend.text = element_text(size=13), legend.title = element_blank(), legend.position = "right") + 
  theme(axis.title=element_text(size=13), axis.text = element_text(size=12)) +
  theme(panel.spacing = unit(0.5, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 11))#, strip.background = element_blank(), strip.placement = "outside")
#annotate("text", label = "2067", x = 2067, y = 60, ymax=80, size = 4, colour = "red") #repeat the y scale for each facet
tag_facet(p0, x = -Inf, y = -Inf, vjust =14, hjust = -0.5, open = "", close = "", 
          fontface = 1, size = 4,  angle = 90, tag_pool = c("2067","2067"))

ggsave(paste0(wd1,"/figures/supplementary/Fig1d.GlobalTrend_Moisture-label2.svg"), device ="svg", units = "in", width = 7.5, height = 2.5, dpi = 600) #horizontal panel

















